function jacobUpdate ()
% Update Jacobian matrix according to matrix "JacBase"
%
global idx_gx idx_gy EXd EXdTJ                % by initial.m
global Gld_cut Bgen_cut Bld_cut               % by initial.m
global busCut Y3fun                           % by Ymatrix.m
global uGall iG_cut signG
global prob nlp                               % nlp.func is modified
global cols JacBase locJx                     % by consJacBase.m

global jpi
if strcmp(nlp.jacobian, 'analytic') == 1
    nlp.func.jacobian = @ jacobianFunc;
else
    jpi = getjpi(nlp.func.constraints, nlp.nvar);
    nlp.func.jacobianstructure = @() jpi.Jc;
    nlp.func.jacobian = @ jacobianAD;
end        

%% Functions
    function Ja = jacobianAD (x)
        [~, Ja] = evalj(nlp.func.constraints, x, [], nlp.ncon, jpi);
    end
    function J = jacobianFunc (x)
        J = JacBase;
        if nargin == 0, return; end
        
        uG = x(nlp.uGidx); uL = x(nlp.uLidx);
        nG = prob.nG; nB = prob.nB;
        uGall(iG_cut) = uG;
        nnode = nlp.nnode; onevec = ones(nnode, 1);     
        % Section <1> <3>
        base3 = (3*nG + 2*prob.ncutG + 4*prob.ncutL) * nnode; 
        idxJ3 = base3-2*nnode+1 : base3;
        for iomg = 1 : nG
            dlt = x(cols(:, iomg)); uEXdi = -EXd(iomg) * (1 - uGall(iomg));
            sindlt = sin(dlt); cosdlt = cos(dlt);
            % Section <1>:  d[F(Domg,...)]/d[delta, Vgx, Vgy]  [NZ = 3*nG * nnode]
            Vgix = x(cols(:, 2*nG+idx_gx(iomg))); Vgiy = x(cols(:, 2*nG+idx_gy(iomg)));
            dVgix = EXdTJ(iomg)*sindlt; dVgiy = -EXdTJ(iomg)*cosdlt;   % d[F]/d[Vgx, Vgy]
            ddlt = -Vgix .*dVgiy + Vgiy .*dVgix;                       % d[F]/d[delta]
            s = 1+ (iomg - 1) * 3*nnode; e = s-1 + 3*nnode; loc = locJx(s : e);
            J(loc) = [ddlt; dVgix; dVgiy];
            % Section <3>:  d[F(V,...)]/d[delta]  [NZ = 2*nG * nnode]
            idxJ3 = idxJ3 + 2*nnode; loc3 = locJx(idxJ3);
            J(loc3) = [uEXdi *cosdlt; uEXdi *sindlt];
        end
        % Section <2>:  d[F(V,...)]/d[V]  [NZ = (2*ncutG + 4*ncutL) * nnode]
        Y3cor = feval(Y3fun, uG, uL);                                  % @ Ymatrix/Y3fun ()
        nzYx = length(Y3cor); idxJ = e-nnode+1 : e;
        for iyx = 1 : nzYx
            idxJ = idxJ + nnode; loc = locJx(idxJ);
            J(loc) = Y3cor(iyx) *onevec; 
        end
        e = idxJ3(end);
        % Section <4>:  [F(V,...)]/d[uG, uL]  [NZ = 2*nparam * nnode]
        idxJ = e-2*nnode+1 : e;
        for ip = 1 : prob.nparam
            Vcix = x(cols(:, 2*nG +busCut(ip))); 
            Vciy = x(cols(:, 2*nG +nB +busCut(ip)));
            idxJ = idxJ + 2*nnode; loc = locJx(idxJ); 
            if ip <= prob.ncutG                                        % [F(V,...)]/d[uG]
                Bi = Bgen_cut(ip); dlt = x(cols(:, ip));
                if signG(ip) == 0
                    J(loc) = [Bi *Vciy ; -Bi *Vcix];
                else
                    J(loc) = [Bi *Vciy + EXd(ip) *sin(dlt); ...
                     -Bi *Vcix - EXd(ip) *cos(dlt)];
                end
            else                                                       % [F(V,...)]/d[uL]
                Gi = Gld_cut(ip -prob.ncutG); Bi = Bld_cut(ip -prob.ncutG); 
                J(loc) = [-Gi *Vcix + Bi *Vciy; -Bi *Vcix - Gi *Vciy];
            end
        end
    end % END of jacobianFunc ()

end